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Abstract 

In  this  paper  we  consider  the  boundary  over  distance  preconditioner  for  radial  basis 
function  interpolation  problems.  We  give  both  theoretical  and  numerical  results  indicating 
that  it  performs  extremely  well. 


1 Introduction 

Let  $ : lZd  — » 1Z,  X = pi, . . . ,x^}  be  a set  of  TV  distinct  points  in  TZd  and  / be  a real 
valued  function  which  we  can  evaluate  at  least  at  the  xfs.  Define 

9-9  = Eili  - xi ) 

where  JV=1  Xjq(x.j)  = 0,  for  all  q G trf 


We  consider  the  problem  of  finding  an  element  s of  + 7 rf  satisfying  the  interpolation 

conditions 


s(xi)  — f(xi),  for  all  xt  € X. 


(1.2) 


Assume  $ is  strictly  conditionally  positive  definite  of  order  2 (SCPD2)  and  X is  uni- 
solvent for  Trf.  Then  there  is  a unique  element  of  S,p^x  + satisfying  the  interpola- 
tion conditions  (1.2).  This  setting  includes  popular  choices  of  the  basic  function  such 
as  the  thin-plate  spline,  $(•)  = | - | 2 log  | - | , and  minus  the  ordinary  multiquadric, 
$(•)  = -\/|  • |2  + c2.  In  this  paper  we  consider  various  ways  of  formulating  the  in- 
terpolation problem,  showing  in  particular  that  a certain  inexpensive  change  of  basis 
can  dramatically  improve  its  conditioning. 

The  usual  way  to  formulate  this  problem  is  in  terms  of  the  functions  {<£(•  — aq)}  and 
some  basis  {po,Pi,  ■ ■ ■ ,Pd}  for  nf.  Then  the  interpolation  conditions  together  with  the 
side  conditions  taking  away  the  extra  degrees  of  freedom  introduced  by  the  polynomial 
part  can  be  written  as 

A\  + Pc  = f and  PT  A = 0,  (1.3) 


where  Ay  = 3>(xj  — Xj), 
4,  5]  that  the  matrix 


Pij  = Pj(xi),  and  / = [/pi),  • • • > /Pjv)]T-  It  is  wel1  known  [3, 


A$  = 


A P 
PT  O ’ 


(1.4) 
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of  this  usual  formulation  is  frequently  badly  conditioned,  even  when  the  number  of  nodes 
is  small.  Indeed  many  authors  have  commented  on  the  numerical  difficulties  that  solving 
this  system  presents  [3,  4,  5] . Results  of  Narcowich  and  Ward  show  that  conditioning  of 
the  system  (1.4)  depends  very  heavily  on  the  geometry  of  the  nodes.  However,  frequently 
in  numerical  analysis  a change  of  basis,  or  other  reformulation,  can  make  a highly  in- 
tractable problem  tractable.  Hence,  our  goal  is  to  find  an  inexpensive  but  highly  effective 
preconditioner  for  RBF  interpolation  systems. 

In  this  paper  we  establish  properties  of  a preconditioning  method  for  the  RBF  inter- 
polation equations  which  was  first  presented  in  Sibson  and  Stone  [5].  In  the  following 
section  we  give  a detailed  account  of  the  preconditioning  method.  In  Section  3 we  prove 
that  the  construction  produces  an  Nx  (N—3)  matrix  Q whose  columns  are  orthogonal  to 
P,  and  which  is  of  full  rank  whenever  the  nodes  X are  unisolvent  for  -k\  . Finally,  Section 
4 contains  numerical  results  for  different  SCPD2  basic  functions  over  a range  of  data 
sets  and  scales.  These  numerical  results  show  that  using  this  inexpensive  0(N  log  N) 
flop  preconditioner  and  variants  of  it,  dramatically  improves  the  conditioning  of  RBF 
interpolation  problems.  See  Figure  1 below. 


(a)  Multiquadric  basic  function.  (b)  Thin-plate  spline  basic  function. 

Fig.  1.  Sorted  2-norm  condition  numbers  of  the  unpreconditioned  matrices,  Aq,.  (top) 
and  of  the  preconditioned  matrices,  S,  (bottom)  for  fifty  thousand  random  data  sets  of 
size  one  hundred. 


2 A preconditioning  method 

A general  approach  to  preconditioning  interpolation  problems  with  SCPD2  basic  func- 
tions in  1Z2  [1,  5]  is  to  choose  Q as  any  N x (N  — 3)  matrix  whose  columns  are  orthogonal 
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to  P and  has  rank  N — 3.  Letting  A = Q/i  and  premultiplying  (1.3)  by  QT  gives  the  new 
system  to  be  solved  for  p,  or  equivalently  A, 

Bp  = QT  f where  B — QTAQ.  (2.1) 

The  three  polynomial  coefficients  can  then  be  found  by  a small  subsidiary  calculation. 

In  this  section  we  present  the  boundary  over  distance  method  of  Sibson  and  Stone  [5] 
for  constructing  the  matrix  Q.  We  will  prove  in  the  subsequent  section  that  Q has  full 
rank  and  is  orthogonal  to  P for  any  set  of  distinct  nodes  X = {xi, . . . ,xn}  C P2, 
which  are  unisolvent  for  n2.  These  properties  of  Q are  well  known  (see  e.g.  [1,  5])  to 
imply  that  the  matrix  of  the  preconditioned  system  B = QT AQ  is  positive  definite.  The 
construction  of  Q is  appealing  in  that  for  “interior”  points  Xj  of  X it  is  local.  That  is, 
for  such  points  the  entries  in  the  j-th  column  of  Q depend  only  on  the  geometry  of  the 
nodes  near  Xj  and  not  on  any  properties  of  nodes  far  away. 

Choose  a closed  bounded  convex  polygonal  region  W of  P2  such  that  X C W . 
Suppose  without  loss  of  generality  that  {xn-2,xn~i,xn}  is  unisolvent  for  it2.  We  will 
refer  to  these  points  as  special  points.  They  are  generally  chosen  so  that  they  are  well 
spread  throughout  W.  In  our  experience,  and  that  of  Sibson  and  Stone,  for  typical  data 
sets  the  choice  of  special  points  is  not  at  all  critical,  as  long  as  the  triangle  they  define 
has  largish  area.  However,  for  contrived  data  sets,  such  as  all  but  a very  few  points 
on  a straight  line,  the  choice  of  special  points  becomes  important.  In  these  cases  we 
have  observed  that  bad  choices  of  special  points  can  lead  to  large  condition  numbers. 
However,  the  strategy  of  choosing  the  three  special  points  to  maximise  the  area  of  the 
corresponding  triangle  has  always  led  to  small  condition  numbers. 

The  region  W is  divided  into  panels  by  intersecting  a Voronoi  diagram  of  the  points 
of  X with  the  region  W.  We  denote  this  panelling  of  W by 

N 

vw{X)  = \J  ^ 

i= 1 

where  V)  is  the  Voronoi  panel  about  the  ith  centre  and  is  defined  by 


V{  = {x  eW  :\x  - Xi\  <\x  - Xj |,  for  all  1 < j < N with  j £ i}. 

Recall  that  the  locus  of  points  equidistant  from  two  fixed  points  is  the  perpendicular 
bisector  of  the  segment  connecting  the  points.  It  follows  that  each  Voronoi  region  is 
polygonal.  Associated  with  a panel  Vi  are  its  edges.  These  are  a finite  number  of  dis- 
tinct closed  line  segments  of  non-zero  length.  They  are  the  boundaries  between  different 
Voronoi  panels,  or  between  a Voronoi  panel  and  Wc.  The  collection  of  all  edges  of  all 
the  Voronoi  panels  will  be  denoted  by  E. 

Definition  2.1  Two  polygonal  regions  of  V?  will  be  said  to  be  strongly  contiguous  if 
they  have  a common  boundary  of  non-zero  length. 

Definition  2.2  Two  Voronoi  regions  Vi  and  Vj  will  be  said  to  be  C-related  if  there  is 
a sequence 


{Vi,  Vtl  ,vt2,...,  Vtm , Vj),  1 < i,  j,  ii,  • • • , tm < N - 3, 
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in  which  all  adjacent  pairs  are  strongly  contiguous. 

Loosely  speaking  V.  and  Vj  are  C-related  if  they  are  connected  by  a chain  of  strongly 
contiguous  pairs.  C-related  is  an  equivalence  relation  on  the  set  {V}^3  °f  Voronoi 
regions  of  non-special  points.  Therefore  it  breaks  this  set  into  a finite  number  of  nonempty 
equivalence  classes  {Qi  : 1 < l < k}. 

Lemma  2.3  Let  Qe  be  any  of  the  equivalence  classes  above.  Then  there  is  at  least  one 
Voronoi  region  Vi  in  Qe  which  is  strongly  contiguous  to  either  Wc  or  one  of 
{Vn-2,Vn-i,Vn}. 

Proof:  Consider  T = Vi  . 

i.ViSQe 

This  union  is  a closed  bounded  connected  polygonal  set  whose  boundary  can  be  written 
as  the  union  of  some  of  the  line  segments  from  E.  Recall  in  particular  that  all  these 
line  segments  have  non-zero  length.  Pick  one  line  segment  < a,  b > from  the  boundary 
of  T.  Since  it  forms  part  of  the  boundary  of  T on  one  side  of  it  lies  a Voronoi  region 
Vi.  from  Qf  . On  the  other  side  lies  either  Wc  or  another  Voronoi  region  Vj.  In  the  first 
case  the  Lemma  is  proven.  Consider  the  second  case.  If  1 < j < N — 3 then  V is 
strongly  contiguous  to  Vj.  Consequently,  Vj  E Qe-  This  contradicts  < a,  b>  being  on  the 
boundary  of  T.  Hence,  N - 2 < j < N and  the  Lemma  follows.  □ 

We  now  detail  the  construction  of  the  N x(N-  3)  matrix  Q using  boundary  over  dis- 
tance weights.  Note  that  because  most  elements  of  Q are  zero  sparse  storage  of  Q requires 
only  0(N)  memory.  A non-special  point  from  {xt  : 1 < i < V-3}  which  has  Voronoi  tile 
that  is  strongly  contiguous  to  Wc  will  be  called  a Voronoi  external  point.  Define  Ve(X') 
as  the  set  of  indices  of  all  Voronoi  external  points.  All  other  points  are  referred  to  as 
Voronoi  internal  points.  The  corresponding  indices  are  Vi(X)  = {1, . . . , N - 3}  — Ve(X). 

We  first  consider  forming  a column  of  Q for  an  index,  j,  such  that  j E Vj(X).  In 
this  case  the  panel  Vj  shares  non-trivial  edges  only  with  other  Voronoi  panels  and  not 
with  Wc . The  column  is  formed  using  boundary  over  distance  weights,  found  from  the 
Voronoi  diagram.  For  j E Vi(X)  the  boundary  over  distance  weight  r.tj  is 


b(xj,Xj) 
\xi  Xj  | 


for  all  Vi  strongly  contiguous  to  Vj, 


(2.2) 


where  b(xi,Xj)  is  the  length  of  the  boundary  between  Vj  and  Vj.  For  other  values  of 
i ^ j,  rij  is  set  to  zero.  In  order  that  column  j of  Q is  orthogonal  to  constants  the 
diagonal  element  rjj  is  specified  as 


rjj  — y ] rij ■ (2-3) 

Finally,  the  jth  column  of  R is  scaled  by  dividing  by  the  area  of  Vj  to  obtain  the  jth 
column  of  Q.  Note  that  the  column  is  by  construction  diagonally  dominant,  but  not 
strictly  so. 

If  j € Ve{X)  then  Vj  is  stongly  contiguous  to  the  complement  of  W,  Wc.  The 
boundary  segment  corresponds  to  a Voronoi  edge  between  Xj  and  an  artificial  point,  the 
reflection  of  Xj  in  the  boundary  (see  Figure  3 in  [7]).  The  reflected  point,  Xj,  can  be 
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written  as  a linear  combination  of  the  special  points,  i.e., 

Xj  — XnXn  + ^N-l^N-l  + ^N-2^N-2,  (2.4) 

where  Aw  + Xn-i  + A jv-2  = 1.  If  Vj  has  k edges  with  Wc  then  k reflected  points 
{£j,.  ■ ■ ,Xj}  are  required.  Associated  with  each  reflected  point,  Xj,  are  the  coefficients 
{Ajv,  Ajy_1(  A^_2}.  The  boundary  over  distance  weights  for  x°  are  partitioned  amongst 
the  special  points  to  obtain  for  all  j G Ve{X)  and  i ^ j 


» Vl  strongly  contiguous  to  Vj, 


(2.5) 


Of  course,  Vj  could  be  strongly  contiguous  with  a Voronoi  panel  associated  with  a special 
point.  If  this  is  the  case  rij  — -M  ■ Again,  for  other  values  of  i / j. 

is  set  to  zero.  Finally  rjj  is  specified  as  in  (2.3)  and  column  j of  Q is  defined  as 


1 ij 


column  j of  R scaled  by  dividing  by  the  area  of  Vj . 
Partition  Q as 

q='e~ 


(2.6) 


where  E is  (N  - 3)  x (TV  - 3).  Thus  E results  from  interactions  between  non-special 
points,  and  F those  between  special  and  non-special  points.  Note  in  the  construction 
above  that  for  1 < i,  j < N — 3,  ev  is  non- zero  if  and  only  if  Vt  is  strongly  contiguous 
to  Vj.  Furthermore,  note  that  E is  necessarily  column  diagonally  dominant,  with  strict 
dominance  in  column  j whenever  Vj  is  strongly  contiguous  to  the  Voronoi  region  of  a 
special  point,  or  to  Wc. 

Relabelling  if  necessary  we  can  assume  the  indices  of  the  Voronoi  regions  in  each  of  the 
equivalence  classes  Gi  form  a contiguous  subset  of  {1, . . . , N - 3}.  Similarly,  we  can  also 
assume  that  the  indices  corresponding  to  any  Gi  precede  those  corresponding  to  Gi+i- 
Furthermore,  by  construction  none  of  the  regions  in  Gi  is  strongly  contiguous  with 

a region  in  Gj ■ Thus,  corresponding  entries  in  the  matrix  E constructed  using  boundary 
over  distance  weights  and  artificial  points  are  zero.  That  is  E is  block  diagonal  with 
the  square  matrix  Eu  on  the  main  diagonal  corresponding  to  the  equivalence  class  of 
Voronoi  regions  Gi-  More  precisely,  Q will  have  form 


' En 

o 

...  o' 

O 

E22 

...  0 

O 

o 

■ ■ ■ Ekk 

. Fi 

f2 

■■■  Fk  _ 

3 Properties  of  the  matrix  Q 

In  this  section  we  establish  the  fundamental  properties  of  the  matrix  Q of  (2.7).  Namely 
that  it  is  of  full  rank  and  that  its  columns  are  orthogonal  to  those  of  P. 
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Definition  3.1  For  m > 2,  an  my.  m matrix  K is  irreducible  if  there  does  not  exist 
an  m x m permutation  matrix  P such  that 

PKPr  = [ l, 

0 M22 

where  Mn  is  r x r,  M22  is  (m  — r)  x (to  — r),  and  1 < r < m. 

The  following  result  is  well  known,  see  for  example  Varga  [6] . 

Theorem  3.2  Suppose  the  square  matrix  K is  irreducible  and  row  (column)  diagonally 
dominant  with  strict  row  (column)  diagonal  dominance  in  at  least  one  row  (column). 

Then  K is  invertible. 

Lemma  3.3  Let  X be  a finite  set  of  distinct  points  unisolvent  for  nf.  Let  Eu  be  one  of 
the  square  blocks  from  the  diagonal  of  Q constructed  in  the  previous  section.  Then  Eu 
is  invertible. 

Proof:  Prom  the  construction  Eu  is  column  diagonally  dominant.  Furthermore,  by 
Lemma  2.3  the  diagonal  dominance  is  strict  for  at  least  one  column  of  Eu.  Prom  the 
definition  of  the  equivalence  relation  C-related  there  is  a chain  of  strongly  contiguous 
pairs  of  Voronoi  regions,  connecting  any  two  Voronoi  regions  in  Qu  This  implies  the  cor- 
responding entries  in  Eu  are  non-zero  and  hence  from  [6]  Theorem  1.6  Eu  is  irreducible. 

It  follows  from  Theorem  3.2  that  Eu  is  invertible.  □ 

Theorem  3.4  The  matrix  Q described  in  Section  2 is  orthogonal  to  P i.e.  QTP  = O. 
Proof:  Omitted,  see  [2]  and  [7].  □ 

Theorem  3.5  Let  X be  a set  of  distinct  points  unisolvent  for  ir'f.  Let  Q be  formed  by 
the  construction  in  Section  2 and  Ay  = — Xj)  where  $ is  strictly  conditionally 

positive  definite  of  order  2.  Then  B = QT AQ  is  positive  definite. 

Proof:  Prom  Lemma  3.3  each  of  the  matrices  Eu  occurring  in  the  block  partitioning  of 
Q given  in  Equation  (2.7)  is  invertible.  Hence  Q has  full  rank.  Also  from  Theorem  3.4  the 
columns  of  Q are  orthogonal  to  the  columns  of  P.  Let  p be  any  non-zero  vector  in  1ZN~3, 
and  define  A = Qp.  Then  A ^ 0,  PT A = PTQp  = 0,  and  pT  Bp  = pTQT AQp  = ATAA. 
Hence,  by  the  definition  of  strictly  conditionally  positive  definite,  pTBp  > 0 whenever 
/i^O  and  B is  symmetric  positive  definite.  □ 

Theorem  3.6  Let  $ be  strictly  conditionally  positive  definite  of  order  2 and  such  that 
$(hx,hy)  = h7<f>(:r,y)  +Ph(  x — y),h>0  with  Ph  € 7rf.  The  preconditioned  matrix  Bh, 
which  corresponds  to  preconditioning  on  the  point  set  hX,  is  a homogeneous  function 
of  scale.  Thus  its  condition  number  and  the  relative  clustering  of  its  eigenvalues  are  the 
same  over  all  scales. 

Proof:  Omitted,  see  [7]. . □ 

Theorem  3.6  applies  in  particular  to  the  usual  thin-plate  spline,  $(•)  = | • |2  log  | • |,  in 
7 72. 

The  extended  version  of  this  paper  [7]  contains  a proof  that  the  elements  By  decay 
like  \xi  — Xj \~K  when  | Xi  — xj\  is  large.  For  the  multiquadric  k is  three  and  for  the 
thin-plate  spline  k is  two. 
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Definition  3.7  The  preconditioned  matrix  S is  obtained  from  B by  pre-multiplying  and 
post-multiplying  B by  the  diagonal  matrix  D with  ii  entry  1 / y/bfi. 

4 Numerical  results 

In  this  section  we  present  numerical  results  for  the  thin-plate  spline  and  multiquadric 
basic  functions.  In  the  following  tables  the  matrix  Aq,  is  defined  in  (1.4),  B in  (2.1), 

5 in  Definition  3.7  and  the  homogeneous  matrix,  C,  is  presented  in  [1].  In  Table  1 we 
show  2-norm  condition  numbers  of  matrices  for  the  various  preconditioning  techniques 
over  seven  different  scales.  It  is  clear  that  the  algorithm  in  Section  2 gives  a matrix 
which  dramatically  improves  the  conditioning  of  the  interpolation  problem.  In  one  case 
by  a factor  of  1014!  Tables  2 and  3 contain  condition  numbers  of  the  matrices  resulting 
from  applying  the  preconditioning  techniques  of  this  paper  for  the  thin-plate  spline  and 
multiquadric  basic  functions.  For  N < 3200,  the  entries  in  the  tables  are  the  maximum 
over  one  hundred  random  point  sets  of  size  N.  For  N = 3200,  the  tables  contain  the 
maximum  over  twenty  random  point  sets  of  size  3200.  In  all  cases  the  preconditioning 
results  in  a smaller  condition  number.  For  these  basic  functions  the  maximum  observed 
condition  number  of  the  scaled  preconditioned  matrix,  S,  grows  very  slowly  with  N. 
Certainly  there  is  no  numerical  evidence  of  power  growth  with  N. 


Scale  parameter 
a 

Conventional 
matrix  A $ 

Homogeneous 
matrix  C 

Preconditioned 
matrix  B 

Scaled 
matrix  S 

0.001 

1.531(11) 

1.534(5) 

4.905(1) 

2.405(1) 

0.01 

1.544(9) 

1.534(5) 

4.905(1) 

2.405(1) 

0.1 

1.597(7) 

1.534(5) 

4.905(1) 

2.405(1) 

1 

3.107(5) 

1.534(5) 

4.905(1) 

2.405(1) 

10 

1.915(6) 

1.534(5) 

4.905(1) 

2.405(1) 

100 

1.271(11) 

1.534(5) 

4.905(1) 

2.405(1) 

1000 

4.006(15) 

1.534(5) 

4.905(1) 

2.405(1) 

TAB.  1.  Condition  numbers  for  one  hundred  points  in  [0,  a]2  and  the  thin-plate 
spline.  The  point  set  for  scale  a is  Xa  = aX\. 


In  an  attempt  to  rule  out  the  possibility  that  our  numerical  results  were  flukes  due  to 
the  small  number  of  100  experiments  we  also  conducted  50,000  trials  with  random  data 
sets  of  size  100.  The  results  of  these  trials  are  shown  in  Figure  1.  The  maximum  condition 
number,  over  all  trials  with  the  thin-plate  spline,  for  the  matrix  A#  was  1.2465(9),  for 
matrix  C,  1.5750(9)  and  for  matrix  S,  1.8066(2).  In  our  experiments  the  matrix  S is 
always  well  conditioned.  This  held  even  for  geometries  of  centres  for  which  the  matrix 
A(f,  is  very  badly  conditioned. 

To  test  further  the  behaviour  of  S for  “bad”  configurations  of  points  a similar  exper- 
iment was  run  with  one  thousand  trials  of  one  hundred  points  almost  on  a circle.  The 
maximum  condition  numbers  of  the  A matrix,  C matrix  and  S matrix  were  1.2885(9), 
7.2692(8)  and  6.6005(2)  respectively  over  1000  trials.  Even  though  the  Voronoi  regions 


Preconditioning  RBF  interpolation 


Number  of 
data  points 

Conventional 
matrix  A# 

Homogeneous 
matrix  C 

Preconditioned 
matrix  B 

Scaled 
matrix  S 

200 

6.555(7) 

3.068(7) 

1.617(3) 

6.028(1) 

400 

5.675(8) 

3.397(8) 

1.945(3) 

8.946(1) 

800 

1.960(10) 

1.348(10) 

2.034(3) 

9.775(1) 

1600 

1.092(10) 

8.413(9) 

8.099(3) 

1.258(2) 

3200 

4.997(10) 

3.783(10) 

1.261(4) 

1.569(2) 

Tab.  2.  Maximum  condition  numbers  encountered  over  a sample  of  100  random 
point  sets  of  size  N in  [0,  l]2  with  the  thin-plate  spline. 


Number  of 
data  points 

Conventional 
matrix  A$ 

Preconditioned 
matrix  B 

Scaled 
matrix  S 

200 

2.014(8) 

1.532(2) 

4.224(1) 

400 

2.045(10) 

5.932(2) 

7.669(1) 

800 

6.641(10) 

4.559(2) 

5.826(1)  ‘ 

1600 

1.554(10) 

7.025(2) 

5.601(1) 

3200 

2.477(11) 

9.362(2) 

6.280(1) 

Tab.  3.  Maximum  condition  numbers  encountered  over  a sample  of  100  random 

point  sets  of  size  N in  [0,  l]2  with  the  multiquadric  function,  parameter  c = 1/y/N. 

are  long  and  thin  the  matrix  S is  still  well  conditioned! 
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